clear; clc; close all;

% 参数设置
q = 1.3; % 分数阶
P = -2; Q = 2;
x0 = 2; y0 = -2;
N_total = 300; % 轨迹长度
N_transient = 200; % 丢弃暂态

a_vals = linspace(-5, 5, 100);
b_vals = linspace(-5, 5, 100);

FraAppEn_map = zeros(length(b_vals), length(a_vals));

parfor ia = 1:length(a_vals)
    a = a_vals(ia);
    temp_row = zeros(1, length(b_vals));
    for ib = 1:length(b_vals)
        b = b_vals(ib);
        [~, y] = SCLMM(q, a, b, P, Q, x0, y0, N_total);
        ydata = y(N_transient+1:end);
        
        m = 2;
        r = 0.2 * std(ydata);
        temp_row(ib) = fractionalApproxEntropy(ydata, m, r);
    end
    FraAppEn_map(:, ia) = temp_row';
end

% 绘图
figure;
imagesc(a_vals, b_vals, FraAppEn_map);
set(gca, 'YDir', 'normal');
colormap(jet);
colorbar;
xlabel('a');
ylabel('b');
title(sprintf('Fig.13(a) 2D 彩色 FraAppEn 分布图 (q=%.2f)', q));
axis square;

%% 分数阶近似熵函数示例
function apEn = fractionalApproxEntropy(data, m, r)
    
    N = length(data);
    phi = zeros(1,2);

    for k = m:m+1
        C = zeros(1,N-k+1);
        for i = 1:N-k+1
            temp = data(i:i+k-1);
            count = 0;
            for j = 1:N-k+1
                if max(abs(temp - data(j:j+k-1))) <= r
                    count = count + 1;
                end
            end
            C(i) = count / (N-k+1);
        end
        phi(k - m + 1) = (1/(N-k+1)) * sum(log(C));
    end

    apEn = abs(phi(1) - phi(2));
end

